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ABSTRACT 


Subsonic and transonic steady and unsteady flowfields over airfoils are inves- 
tigated with the numerical solution of the governing equations. This study aims to 
enhance the performance of rotary wing and fixed wing aircraft by better understand- 
ing and by taking advantage of unsteady phenomena such as dynamic lift. In the past 
few years many advances have been made in algorithm development for the numerical 
solution of the Euler and the Navier Stokes equations. In this study, these new zonal 
techniques are applied. A zonal approach is more computationally efficient in solving 
the governing equations than previous approaches, and has certain advantages over 
the standard single moving grid approach. The zonal grid consists of two grids, one 
being the inner grid which is fixed to the airfoil, and the other being the outer grid 
which extends to the far field or to a specified outer boundary. The inner grid is 
allowed to rotate with the body, while the outer grid remains fixed. The thin-layer 
Navier-Stokes equations are solved for the inner grid, while the Euler equations are 
solved for the outer grid. Communication between the two grids is accomplished by 
interpolating the flow quantities at the zonal interface. Solutions are obtained for 
flows at fixed angles of incidence, and for unsteady flows over pitching and oscillating 
airfoils. The computed results are in good agreement with available experimental 


data. 
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I. INTRODUCTION 


A. BACKGROUND 

Investigation of steady and unsteady flowfields over airfoils is an active area 
of numerical and experimental research. Unsteady pitch-up motion of airfoils alters 
significantly the aerodynamic characteristics of lifting surfaces. Pitch-up motion of 
an airfoil produces lift augmentation and delay of stall to higher angles of incidence 
compared to airfoils held at a steady incidence. Understanding the mechanisms that 
cause the dynamic lift development is a subject of interest in both theoretical and 
applied research. 

A comprehensive explanation of the dynamic lift phenomenon is given by Tyler 
and Leishman [Ref. 1]. To briefly summarize, unsteady airfoil motion tends to 
maintain high lift to higher angles—of-attack, because the unsteady flow causes a 
time delay in the build-up of the lift force and the adverse pressure gradient. The 
unsteady motion also gives the airfoil a virtual camber which decreases the leading 
edge pressure and pressure gradients. Depending on the flow conditions, as the airfoil 
begins to stall, a vortex forms at the leading edge which grows with time which 1s 
eventually shed downstream in the wake. This vortex shedding phenomenon alters the 
chordwise pressure distribution on the upper surface of the airfoil resulting in higher 
maximum lift coefficients. Sometimes. however, only a recirculating flow is observed 
at the airfoil trailing edge. As the pitch-up motion progresses the recirculating region 
extends upstream towards the leading edge and eventually a vortex is formed at the 


trailing edge. 


Rotary wing and fixed wing aircraft designers can enhance the performance of 
the aircraft by taking advantage of the dynamic lift concept. However , the resulting 
undesirable pitching moment variations have deprived helicopters and airplanes of 
the benefits of dynamic lift [Ref. 2]. Carr [Ref. 3] gives a comprehensive review of 
the progress that has been made in the study of the dynamic stall phenomenon. 

Experimental work on steady transonic flow on airfoils was conducted by McDe- 
vitt and Okuno [Ref. 4]. These data cover a wide range of flow conditions and can 
be used for steady code validation. One of the pioneering experimental works on un- 
steady airfoil flows was conducted by McCroskey et al [Ref. 5]. In their experiments, 
unsteady dynamic effects for two-dimensional airfoil flows were studied for the first 
time in conjunction with flows over helicopter blades. They showed the effects of 
unsteady lift and pitching moment on the retreating rotor blades. Their research 
also showed that unsteady motion parameters such as reduced frequency, appeared 
to be more important than airfoil shape in determining the dynamic stall airloads. 
Benchmark experimental data of oscillating and pitching airfoils was also collected by 
Landon [Ref. 6]. The experimental data of references [5] and [6] gives enough quali- 
tative information for unsteady code validation. Recent experiments with oscillating 
airfoils, performed by Chandrasekhara and Brydges [Ref. 7], have shown definitively 
that at subsonic Mach numbers the unsteady flow in the vicinity of the leading edge 
can reach supersonic speeds and generate shocks. 

Along with all the experimental studies being conducted, a great effort is also 
underway to compute unsteady viscous flows. Developments of numerical methods 
for the Navier-Stokes equations [Refs. 8, 9, 10] during the past few years provide 
new tools for the investigation and prediction of airfoil flows. In this investigation 
the compressible thin layer Navier-Stokes equations are solved using a zonal grid 


approach. The objective is to develop a computationally efficient method to study 


い ウ 


two-dimensional unsteady flows. This will be accomplished by developing a solution 
procedure for two grids, an inner viscous grid around the solid body, and an outer 
grid which is coarser representing the outer flow field. The inner grid is free to rotate 
to any angle of attack. and the outer grid remains stationary. 

The idea of moving meshes and dynamic meshes is not new. The dynamic 
and adaptive grid solution method refers to computational grids which are coupled 
to the physical problem that is being solved.An adaptive solution procedure with 
grid points that continually move during the solution process in order to resolve the 
developing gradients has been shown [Ref. 11]. A dynamic type of mesh has been 
applied by Rumsey and Anderson [Ref. 12] to simulate aileron buzz using the thin 
layer Navier-Stokes equations. 

The idea of overlapped or patched grid schemes was used by Rubbert and Lee 
[Ref. 13] with the limitation that the grid lines at the boundaries were continu- 
ous. Benek et al [Ref. 14] developed the “CHIMERA” approach for two-dimensional 
embedded and overlapping grids. They demonstrated that the grid lines at the over- 
lapped boundaries did not have to be continuous, and that the flow quantities could 
be successfully interpolated. The “CHIMERA” approach was later extended to three 
dimensions by Benek et al [Ref. 15] . Rai [Ref. 16] developed a technique for in- 
dependent zonal grids where the flow variables were interpolated across the zonal 
interface. Chesshire and Henshaw [Ref. 17] developed a methodology for solving the 
steady compressible Navier-Stokes equations using multiple overlapped grids. Their 
methodology was demonstrated by solving the governing equations on a composite 
grid that was comprised of an airfoil grid, a leading-edge flap grid, a trailing-edge 
flap grid and a grid for the outer flowfield. Chyu and Davis [Ref. 18] investigated 
unsteady transonic flow using a moving airfoil grid. They developed stationary com- 


putational grids around the airfoil at its lowest and highest angles of attack, and then 


interpolated a new grid as the airfoil oscillated to intermediate angles—of-attack. Reu 
and Ying [Ref. 19] developed a composite grid approach to study the flow about 
pitching airfoils in a wind tunnel. Their approach consisted of a structured inner grid 
and an unstructured grid for the outer flowfield. 

Unsteady problems have also been solved by oscillating the incoming flow, as 
opposed to moving the grid. However, this method does not work if a solution is 
sought for two or more objects, such as a wing tail combination or a canard wing in 


relative motion to each other. 


B. PURPOSE 

In this work, unsteady compressible flows are investigated using a numerical 
technique which is applied to zonal meshes. The governing equations are solved 
on multiple computational grids, where one of the grids is free to move in unison 
with the solid boundary and the other grid is fixed. The meshes overlap at the 
zonal interface. The scheme that is developed is similar to the “CHIMERA” scheme 
mentioned previously, but we impose the restriction of a circular shape on the zonal 
interface. 

The new approach developed in this study is more computationally efficient 
compared to previous schemes used to study unsteady flows. The zonal grid approach 
avoids the need to regenerate or interpolate entire grids at every angle of attack. The 
zonal interface is in a known position; therefore the entire flowfield does not need to 
be searched. In addition, since the outer grid remains stationary, the metric terms do 
not need to be recomputed at each time step. Another advantage is that the zonal 
interface is circular; therefore the flow variables are interpolated in the circumferential 
direction only. The zonal grid approach also allows for the application of different 


solution methods to the inner and outer grids. For example, viscous solution methods 


near the solid body and an inviscid method on the outer mesh can be implemented. 
The primary objective of this investigation 1s to develop and test the zonal grid 
methodology. The space discretization is based on Osher’s [Ref. 20] upwind method. 
An implicit scheme is used for time integration. An advantage of upwind schemes is 
that they are naturally dissipative and no explicit artificial dissipation is required. 
In order to meet our goal, a procedure for generating zonal grids was devel- 
oped, along with the zonal flow solver. The effect of grid resolution on the accuracy 
of solutions was studied first for steady flows and the solutions were compared to 
experimental data collected by Harris [Ref. 21]. The unsteady results were validated 
using the experiments by Landon. The dependence of the solution on the location of 
the overlapped zonal region was examined. Viscous steady state solutions were com- 
puted and compared to experimental data. Finally, unsteady flows were investigated 
by computing solutions for airfoils in ramp and oscillatory motion. The ramp motion 
was started at zero angle of incidence. and ramped up to 30 degrees. The computed 
pressure coefficients were verified by comparing with available experimental pressure 
distributions. The boundary layers computed for this case were also compared to an 
interactive boundary layer code [Ref. 22]. The oscillatory test case was verified by 


comparing the computed pressure coefficients with experimental data. This approach 


and the computed results are described in the following chapters. 


Il. GOVERNING EQUATIONS 


In order to compute compressible viscous fluid flow around a body, the conti- 
nuity, momentum and energy equations must be solved simultaneously. The vector 
and the conservation-law form of the compressible Reynolds—averaged Navier-Stokes 


equation is presented. A detailed derivation can be found in (Ref. 23] 


A. CONTINUITY EO tits 


The continuity equation expresses the conservation of mass law applied to a 


fluid passing through a control volume fixed in space 


Ор 
— . pV) = 2.1 
ጋው ከከ 0 (2.1) 


here p is the fluid densitv and V is the fluid velocity. Equation 2.1 states that the 
net mass flux through a control volume bounding surface must be equal to the time 
rate of change of the mass inside the control volume. In a two-dimensional Cartesian 


coordinate system this equation reads 
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where u and w are velocity components along the z and z directions, respectively. 


B. MOMENTUM EQUATION 

The momentum equation expresses Newton’s second law as applied to a fluid 
element passing through a control volume fixed in space. The momentum equation 
1S: 


д 
ААН ХАЛАГ ААН (2.3) 


The first term in equation 2.3 represents the time rate of change of momentum 
per unit volume in the control volume. The second term represents the moment flux 
through the bounding surface of the control volume. f is the body force per unit 


volume and IL, 18 the stress tensor given by 


Qu; 2 
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da, 3 995, 
where 2,7,4 = 1,2,3 and ó, is the Kronecker delta. 
By substituting equation 2.4 into equation 2.3 and expanding equation 2.3 for 


a two-dimensional Cartesian coordinate system obtain 


le 
に 


Ji 


p% = pfe- Z+ э; |1н (25% — Э=)| + #2 |n (29 4 


2 ef ተ 22 18 2፳- ፳)] * x ET) 


Equations 2.5 are known as the Navier-Stokes equations for two-dimensional flow. 


C. ENERGY EQUATION 
The energy equation is derived by applying the first law of thermodynamics ( rate 


of change of energy = net heat flux into particle + rate of work done on particle. ) 
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where e is the total energy per unit volume and Q is the heat addition per unit 


volume. 


The above equations can be rewritten in non-dimensionalized vector form as 
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En Ц Ox Oz Rel Oz Oz eM 


where 


pu pw 
pu? +p puw 

ea, 

puw pw + p 

(e + pju (e+ p)w 
0 0 

Trt G. = 

7 2 T, z 


fa 94 


Here. 


fa = UTrr T UTrz + ሙን ST 
r 


94 = uTrz r UTzz ዣ ! а 


MESE 


U 


Re 





where Uoo ıs the free stream speed and L is a reference length The pressure is related 


to u,w,p, and e by 


= ¿o + 107) 


In the above equations, the ratio of the specific heats, y, 1s set equal to 1.4 and 


a* = yp/p is the local speed of sound. 


The density 1s non-dimensionalized by the free stream density, po, the velocities 


are non dimensionalized by the free stream speed of sound, and the total energy is 


non dimensionalized by p.a%,. 


D. TURBULENCE MODEL 


The Navier-Stokes equations can completely model fluid flow, however, in or- 
der to resolve turbulent scales at high Reynolds numbers aud realistic geometries. 
very high grid densities are required. Therefore. present state-of-the-art algorithms 
and computer technology allows direct simulation of turbulent flows for only simple 
geometries and low Revnolds number flows which are of limited practical interest. 
In order to enable computation of turbulent flows for configurations, of practical in- 
terest, turbulence modeling is used. Turbulence models are implemented with the 
time-averaged forms of the Navier-Stokes equations. 

Two widely used averaging procedures are the standard time averaging proce- 
dure for incompressible flow, and the mass-averaged approach for compressible flows 
‘Ref. 24]. The time averaging procedure destroys high frequency information of the 
turbulence, but the unsteady mean flow information is preserved. 

For the incompressible case, the randomly changing flow variables are replaced 
with their averages plus their fluctuations. For example, in the Cartesian coordinate 
system. the u velocity component is represented as ú = ü + и’. where ü is the mean 
velocity and wu’ is the fluctuation about the mean. The governing equations are time 
averaged, and the average of the fluctuation terms 1s set equal to zero. 

In the compressible flow case, the mass-weighted variable of the Favre averaging 


approach is used. In this case, u is represented as u = ù + u” where ù is 


E 
| 


や | |. 


Here the time average of the doubly primed fluctuating quantities is not equal to zero. 
After the substitution is carried out for all of the fluctuating flow variables, the entire 


equation is time averaged. Next, all the time averaged terms that are doubly primed 
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and multiplied by densitv are defined to be zero. For example 


The equations of mean motion resulting from the time averaging procedure have 
more unknowns than equations. This constitutes the closure problem of turbulence. 
In order to close these equations a turbulence model is used. Several models have 
been proposed, in this study. 

The Baldwin-Lomax (B — L) turbulence model |Ref. 25] was used. This model 
is a two-layer eddy viscosity model which simulates the effect of turbulence in terms 
of the eddy viscosity coefficient u;. The term p in the stress terms is replaced by 
и-+ ui, and the u/ P, in the heat flux terms is replaced with L/P, +pu,/ P. . The B— L 
turbulence model is similar to the Cebeci-Smith turbulence model [Ref. 26]. but it 
bypasses the need for finding the edge of the boundary layer by using vorticity instead 
of the boundary layer thickness. This model is adequate for flows which have mild 
pressure gradients, but it is not very suitable for highly separated flows. A complete 
description of the model is given in reference [25]. The basic equations of the model 
tollow. 

In the inner layer, the eddy viscosity is assumed to be proportional to the mixing 
length and vorticity, and in the outer layer it uses an exponentially decaying formula. 
The inner eddy viscosity is computed up to the point where it is equal to the outer 


eddy viscosity as shown below. 


Ub 
254 
— 


THES (Ht inner y > Ucrossover ( 
| (ie eee Ucrossover < y 


where y is the normal distance from the wall and Yerossover taken at its minimum value 
where it equals y. The inner eddy viscosity is given by: 
2 
(27. = pl m 
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+ 
(= à 一 Es 


a 
= Ue 


y PwUrY 


Нэр 


where 


y 
AT is an experimentally determined damping constant, « is the Von Karman constant 
The outer eddy viscosity 1s given by 


Houter = KC opp F (y wake (y)KLEB- 


F(y)krEB 1s the Klebanoff intermittency factor given by 


CKLEBY | | 


Umax 


I 


F(y)kLEB ^ |1-- 5.5 | 





« and C, are constants. For boundary layers 
Il каке = US LS 
For wakes and separated boundary layers 
UDIE 
A E Co が max コー 
max 


The quantity ymax is the value of y determined for the maximum value of Fmax and 


Fmax 15 determined by 
F(y) = vll [1 ep rl. 


E. TRANSFORMATION TO GENERALIZED COORDINATES 
In order to use an unweighted differencing scheme that facilitate the numerical 


implementation on body fitted coordinate systems suitable for complex configura- 


tions, the equations are transformed to a generalized coordinate system using the 
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following transformations 
a mS (2.8) 


The above equations transform the governing equations from the physical domain to 
a body fitted coordinate system. The transformation is carried out by using the chain 


rule of partial differentiation 


ው. P RN 
Ox cM roc DE Зу 


(2.9) 


For example, the continuity equation would be transformed in the following manner. 


dpu Шан дро дро — gonan T Я 
dr тэг дЕ + сс ፀር ' 9: = E us (235 


Ex, €, (y and С, аге known as the metrics. The metric terms can determined in the 


following manner. First write the differential expressions for £ and Ç 
de = £ dz + Са ma E (2718 


In matrix form the expressions become 


as . Ss | 2 (2.12) 
ас CA dz 
Next we write the differential expression for z and z 
dx = r d£ + rcdG (2187) 


dz 一 z d£ _ z dÇ 


In matrix form they become 


a c dé л 
ИШКЕН ግ 


By comparing equations 2.12 and 2.14 we can write 


i x] -[% 3 
Z =ረ e € 


m 


and solving for the metrics we get 


2 = "s 
fF ==Jx | 
7 : (2.16) 
Cr = =J ze 
Gz 一 y SE 
where 
l | 
po unt 
Te Tc 








The mapping from (r, z) to (€,¢) is one to one if the Jacobian, J. is non singular. 


F. THIN LAYER NAVIER STOKES EQUATIONS 


In order to accurately calculate the normal gradients in the boundary layer. it 
is necessary to make the normal grid spacing very fine close to the solid surfaces. n 
the streamwise direction the flow gradients are not large and no fine grid spacing is 
required. As a result the grid cells near the body have a very high aspect ratio. With 
this type of grid, existing gradients in the streamwise direction are not fully resolved. 
In order to facilitate the numerical implementation of the thin layer approximation 
is employed by retaining only the viscous derivatives normal to the body. The thin 
layer Navier-Stokes equations transformed into a generalized curvilinear system Is 


Esctor form.are 


0Q OF óG 1 В 


[DEN NE Бо 
where 
p 
~ 1 | pu 
e 
pU 


(Ep 
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pW 
ë = | puW - Gp 
E pwW t Gp 
(eap UD 
The viscous flux terni 1s transiermeasas 
0 


g _ Û | amu + (1/3)mó 
J| Om o 
umma  (u/3)mom4 


where 


== E + CG 
Hi = 2-1 + Ny, 


-9 (42-42) 1 242 


e A 
mamoa Tomas I por =m 80 


11218552 21 ин са. 


U and W are the contravariant velocity components. These velocity components 


are normal and parallel to the constant € and (С surfaces, respectively. 
С = E + Gud Ew 


Wee €. + CABE? (0 
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III. SOLUTION METHODS 


A. ZONAL GRID GENERATION 

The zonal grids were generated using a software package called GRAPE2D. de- 
veloped at NASA Ames Research Center by Reese L. Sorenson [Ref. 27]. GRAPE2D 
generates field grids by solving the following set of elliptic Poisson equations: 


Баб 
SIT 


£, = P 
Maz — zz = Q 


(3.1) 


GRAPE2D was used to generate the inner and outer grid separately. In gener- 
ating the inner grid, the number of circumferential and radial grid points, the spacing 
at the body surface, and the radius and shape of the outer boundary were specified. 
Once the inner grid was generated, the outer boundary of the inner grid was used as 
the inner boundary for the outer grid. 

Due to the placement of the reference axis at the quarter chord of the airfoil. 
the grid spacing in front of the airfoil is larger than aft of the airfoil. This causes a 
problem in specifying an initial grid spacing for the outer grid. In order to obtain 
the smoothest overall transition from the inner grid to the outer grid, the starting 
spacing for the outer grid was obtained by averaging the minimum and maximum 
spacing at the outer boundary of the inner grid. The outer boundary of the outer 
grid was prescribed as a rectangular region, usually about six chord lengths away 
from the body. 

The perfectly circular zonal boundary was accomplished by reading the over- 
lapped regions r and z coordinates and finding an average radius which was measured 
from the center of the grid, regardless of the location of the inner body. Next, an 


arbitrary grid point was chosen to use as a reference point. Then an angle theta, 
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0, was calculated which represented the angle between the reference grid point and 
the grid point under consideration. Then having the average radius and the angle 0 
for the grid points in the overlapped region, a new set of x and y coordinates was 


calculated using the following equations: 
тав гс0:09 750 10 (3.2) 


The overlapping of the two grids was performed by adding the next to last layer 
of grid points from the inner grid to the outer grid, and adding the second layer of 
grid points from the outer grid to the outer boundary of the inner grid. This resulted 
in an overlapped region of three grid points or two grid cells in the radial direction. 


The grid was also overlapped by a grid point in the circumferential direction. 


B. NUMERICAL SCHEME 


The numerical integration is performed using an upwind-biased, factorized, 


iterative, implicit numerical scheme [Ref. 20] given by 


La ARAS + МА) х 
7 + he (VB, + Ai Bir = Вет, М.)| x рше = (3.3) 
zo Ж a t (Вук Fi aa) 


~ 


HG a = ET = RS = S k-i/2)] 
In equation 3.3. he = Ar/A€, etc., ÀÁ* = (9F/0Q), etc., are the flux Jacobian 
matrices, and A, V, and ó are the forward, backward and central difference operators. 
respectively. The quantities 20117 Са and Sikri f2 are numerical fluxes. The 
inviscid fluxes F and G are evaluated using Osher’s upwinding scheme [Ref. 19]. The 


numerical fluxes for a third-order accurate scheme are given by 


^ 


ት 12 ES p 2 5 arem + 2A FH jaa ー 


~ 


6 or m 十 “一 250 = Fiir нь E 
a LAF (Que Qu) t 2AP* (Qus Qia)] — 
z[AF-(Quu 077 司 
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Where F is the first-order numerical flux for the explicit Osher’s scheme given by 
z ] | Qi41 + ” 
Fina = 5 [Fiat Erna [ (FP - Fp} aq 
- Qi 


pere [^ = hr - eae E = (E). and AF are the corrections to obtain a 
higher order accuracy. The Osher scheme evaluates the flux assuming a shock tube 
solution where Ё, IS piecewise continuous and yields good predictions of the flux. 
especially at supersonic Mach numbers. For the linearization of the left-hand side of 
Eq. 3.3, the flux Jacobian matrices A, B are evaluated by the Steger- Warming [Ref. 
28| flux-vector splitting. 

lime accuracy of the implicit numerical solution is obtained by performing 
Newton iteration to convergence for each time step. The approximation of Q”*! 
at each subiteration is the quantity QP. When p > 2, during a given subiteration, 
OP — Q^*!. but when p = 1 and no subiteration is performed, then Q? = 0", апа 
QP?*! — Q^*!. By subiterating to convergence, linearization and factorization errors 
are minimized, because the left-hand side of Eq. 3.3, can be driven to zero at eacli 
Ane step. 

The linearızation errors are eliminated by subiteration to convergence. Tvp- 
ıcally, two to three subiterations are sufficient to drop the residuals two orders of 
magnitude during the Newton iteration process. 

High order accurate shock-capturing schemes have some limitations; they may 
select a nonphysical solution, they produce spurious oscillations and they may de- 
velop a nonlinear instability in nonsmooth and discontinuous flow regions. More ap- 
propriate high-order shock-capturing schemes, suitable for the computation of weak 
solutions are the TVD schemes, described in detail in [Ref. 29, 30]. In the present 
study, the Osher-Chakravarthy [Ref. 31] TVD scheme ts used. This TVD scheme has 


flux limiters which impose constraints on the gradients of the fluxes. The flux-limited 


ji 


values \f* are computed from the unlimited fluxes Af* as follows 


Ши = minmod ЗУЛ k> IAS ET, a 


Nf asa 2 minmod ana س‎ a 
дэс ЭГ 
ХА, ууу, лс тилтой 1: Ey fiu NT 2 | 
at : 
Af рро = minmod | \ ft TT MM J 
where the minmod operator is defined by 
minmod|r, y] 2 sign(z) x maz[0, min(|z], ysign(x)]] (3209 


The viscous fluxes S; +1/2 are computed with central differences as follows 


a te 1 [2 m SO ES KU Ci k+1/2| 1200) 
Ї 

бо = 5(Qi4 + Qik-ı) (3.8) 

(Qc)ik+1/2 = Qik = Qik (3.9) 


The experimental Reynolds numbers based on the chord length for the test cases 
examined are in the range Re, zz 3.0 x 106 — 5.0 x 105, and it is expected that the flow 
is mostly turbulent. Transitional flow 1s expected to have a small effect at regions very 
close to the leading edge. Present knowledge about boundary laver transition does 
not enable computation and modeling of the transition regime. Therefore, only fully 
turbulent solutions were computed. In the present work, the widely used two-layer 
Baldwin-Lomax turbulence model was used. The effectiveness of other turbulence 
models, such as the Johnson-King model [Ref. 32] and the RNG based algebraic 
model [Ref. 33] for steady and unsteady flows, was investigated in references [34] and 


[35]. 


C. BOUNDARY CONDITIONS 
The solutions on the two grids are computed separately, with the inner and 


outer solutions communicating through the zonal interface boundary. The inner grid 
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surrounds the airfoil and includes the boundary layer region and the wake for viscous 
solutions. Inviscid solutions are obtained by applying the flow tangency slip condition 
where the normal contravariant velocity component, HW’, is set equal to zero on the 
surface. For viscous solutions, the nonslip condition is applied for the velocities on the 
airfoil surface. In both cases the density and pressure are obtained from the interior 
grid points by simple extrapolation. For unsteady solutions, the surface velocity is 


set equal to the airfoil speed obtained by the prescribed airfoil motion as follows 


| : l . 
= j Gi: ee) = 7 (816: 214209 


Unsteady solutions for pitching and oscillating airfoils are obtained by rotating the 
inner grid only. Therefore, only the metrics of the Inner grid must be reevaluated for 
each time step. 

At the inner zonal interface, the flow variables are obtained from the interior of 
the outer grid solution. Similarly, the inner zonal boundary of the outer grid obtains 
boundary information from the interior of the inner grid. The inner and outer grid 
radial lines are not aligned, in general. The relative location of the two grids with 
respect to the inertial reference frame as the inner grid rotates is computed. These 
distances between neighboring points at the zonal interface are used tor a weighted 
averaging of the conservative variables. 

All flows were computed for subsonic free-stream speeds. At the subsonic inflow 
and outflow boundaries of the outer grid, the flow variables were reevaluated using 
zero-order Riemann invariant extrapolation. At the inflow boundary, there are three 
incoming and one outgoing characteristics. Therefore, three variables, the density p, 
the normal velocity w, and the pressure p, are specified and the fourth variable, the 
axial velocity u is extrapolated from the interior. The inflow boundary conditions are 


given by 
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а= e , — Ву) 
i eee ie 


Ср = WY 


p =( 


ማሣ 


here RT. RZ are the incoming and outgoing Riemann invariants given by 
+ Е | E 
Hl =e, 2ас2/(т-үЇ эт iM 
\t the outflow boundary there are one incoming and three outgoing characteristics. 
|! herefore only one quantity, the pressure, is specified. while the others are extrap- 
ulated from the Interior. For the density and normal velocity. simple first-order 
extrapolation is used, and the axial outflow velocity is obtained from the zero-orde: 


outgoing Riemann invariant. The outflow boundary conditions are given by 


Ara 
uy = RI — 2a1/(7y — 1) 
E ол (93:11 
u = 009 


Dee) 


were solved. For the unsteady flow solutions, the outer grid remained stationary and 


the metrics were not reevaluated at each time step. 


IV. RESULTS AND DISCUSSION 


The validity of the zonal grid approach was first investigated for inviscid solu- 
tions. An advantage of the present approach is that different grid densities may be 
used for the inner and outer grids. However, the accuracy and the conservative char- 
acter of the solution for different locations of the zonal interface and grid densities 
must be assessed. 

First, the accuracy of the computed results for different inner and outer grid 
densities was evaluated. The effect of the location of the zonal interface relative to the 
airfoil on the accuracy of the solution was also investigated. Then viscous solutions at 
fixed angles of incidence, up to approximately the static stall angle, were computed. 
Finally, unsteady flow responses to a ramp motion at subsonic free-stream speed of 


M.. = 0.3 and for an oscillation at a free-stream speed of M. = 0.6 were computed. 


A. STEADY STATE SOLUTIONS 
1. Inviscid Test Cases 

Preliminary test cases were computed using coarse meshes with an inviscid 
solution. A two-block grid consisting of an 81 x 40 point O-type inner grid and an 
31 x 22 point O-type outer grid was used as a baseline grid for the inviscid solutions. 
Table 4.1 gives the inviscid grids that were tested. 

a. Case 1. Baseline Grid: 81x40 Inner and 81x22 Outer 

An inviscid solution using the baseline grid for subsonic flow over a 

NACA-0012 airfoil at \/,, = 0.8,a = —.1° was obtained. The baseline grid is given in 
figure «4.2. The distribution of the computed surface pressure coefficient is compared 


with the measurements of [Ref. 4] in Fig. 4.1. Agreement with the experimental data 


トウ 
LO 


NN GRID) DENSITIES OF THE INVISCID SOLUTIONS COMPUTED 
Meets ACA-0012 AIRFOIL AT Af,, = .8 AND a = —0.1. 


| 1 | 81x40 | 51 х 22 IS CHOC | 
| 2 H x 22 1.5 x chord 
| лога 
OO | 
ГЭ 318 STI Doli | 
6 3:39 Ix | 13 x chord | 


is satisfactory for an Euler solution. It can be seen that with this mesh the strength 




































of the shock is captured, but the location is lagged by 5 percent of the chord. The 


baseline solution is converged at 2000 iterations. 
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Figure 4.1: M5 = 0.8 Computed Solution With Baseline Grid Compared to Exper- 
ment. 
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b. Case 2. Baseline Inner Grid With a 41x22 Outer Grid 
Next the effect of the outer grid density is investigated. The grid for 
Case 2 was generated by starting with the baseline grid and removing every other grid 
point from the outer grid. The resulting grid is given in Fig. 4.4. The Euler solution 
is compared to experiment in Fig. 4.3. Overall the converged solution for this case 
agrees with the baseline solution except it 1s noticed that the shock location for the 


upper and lower surface are slightly farther apart than the baseline grid solution. 
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Figure 4.3: \/,, = 0.8 Computed pressure coefficient solution compared to experi- 
mental data. 
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Figure 4.4: Inner Grid 81 x 40 Outer Grid is 41 x 22 
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c. Case 3. Overlap Boundary Set to 1 Chord Length Away From 

Body 
In this case the effect of the overlap boundary location is studied. The 
grid generated is given in Fig. 4.6. It is seen that the transition from the inner to 
the outer grid is not as smooth as for the cases where the boundary was located 50 
percent farther from the airfoil. At the leading and trailing edges. the overlapped 
boundary is only a half chord length away. In figure 4.5 the converged solution 18 
given. With this grid very little effect on the shock strength and location is seen. The 
pressure coefficient before the shock is slightly underpredicted and after the shock is 


slightly overpredicted. 
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Figure 4.5: M5 = 0.8 Computed Solution Compared to Experiment. 
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Figure 4.6: Inner Grid 81 x 40 Outer Grid is 81 x 22 
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d. Case 4. Overlap Boundary Set to .75 Chord Length Away 

From Airfoil 
In order to further investigate the tendencies observed in Case 3, an- 
other grid 1s developed with the zonal interface closer to the airfoil. The overlapped 
region is only a quarter chord away from the leading and trailing edges of the airfoil 
as seen in Fig. 4.5. The solution obtained is compared to experiment in Fig. 4.7. 
The pressure coefficient is again slightly underpredicted before the shock and slightly 
overpredicted after the shock. It is also observed that the zonal interface location 


relative to the airfoil has little effect on the solution. 
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Figure 4.7: .\/,. = 0.8 Computed Solution Compared to Experiment. 
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Figure 4.8: Inner Grid $1 x 40 Outer Grid is 81 x 22 
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e. Case 5. Grid With Oval Interface 

The solution obtained for cases 3 and 4 showed that the location of 
the zonal interface had very little effect on the computed solutions. The ability of 
the zonal interface to pass flow discontinuities was also studied. It was known from 
the previous cases that the shock location was near the the mid chord point. The 
inner grid for this test case had to be generated so that the overlapped region was 
very close to the upper and lower surface of the airfoil. This was accomplished by 
generating an oval inner grid and an oval zonal interface. This grid is shown in figure 
4.10. The computed surface pressure coefficient distribution is shown in figure 4.9. It 
is in agreement with the experimental data and with the previous computed solutions. 
The computed flow quantities, such as density and pressure, showed that the zonal 
approach used can pass shocks through the zonal interface. Figure 4.11 shows Mach 


contours which smoothly pass through the overlapped boundary. 
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Figure 4.9: Mo = 0.3 Computed Solution Compared to Experiment. 
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Figure 4.10: Inner Grid 81 x 18 Outer Grid is 81 x 19 
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Figure 4.11: Mach Contours: Overlapped Region's Ability to Pass Flow Disconti- 
nuities. 
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f. Case 6. Baseline Grid With Half the Radial Grid Points On 

Outer Grid 
In the previous cases, the effect of circumferential resolution was stud- 
ied. Next. the effect of the radial resolution on the computed solutions is investigated. 
Figure 4.13 shows the grid generated to test these effects. The computed surface 
pressure coefficient distribution (in figure 4.12) was comparable with the solutions 
obtained with denser outer grids. Computations with an even coarser grid, e.g., a 
41 x 11 point grid, predicted the shock location even further downstream due to lack 


of streamwise resolution. 
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Figure 4.12: Af,, = 0.8 Computed Solution Compared to Experiment. 
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2. Viscous Test Case 

The Euler solutions presented in the previous section were used to study 
the effects of different grid densities on the computed solutions. Having confidence 
in the predictions of the code, the next step was to generate a viscous grid with the 
airfoil’s quarter chord point at the center of the inner grid. 

Viscous, subsonic flow solutions were obtained at the following fixed angles 
of attack: 3.27°.4.97°,6.69°78.38°, 9.27°, 10.1227 10.992 And 11.00 The How s 
tions of the measurements reported in [Ref. 21] were used, e.g., Mœ = 0.3. Re = 
4.0 x 109. The spacing of the grid was set to 0.0005 for the first grid point above 
the surface. The quarter chord point of the airfoil was set at the center of rotation 
so the airfoil could be ramped about the quarter chord point. These solutions were 
obtained on a 181 x 56 point viscous inner grid and a 181 x 26 point inviscid outer 
grid shown in figure 4.14. 

Solutions were also computed on a grid with half the streamwise resolution. 
e.g.. à 91 x56 point grid. The computed surface pressure coefficient distributions using 
the 181 x 56 inner grid and 181 x 26 outer grid, are compared to experimental results 
in figures 4.15 through 4.23. 

Solutions for fixed angles of incidence were obtained by two methods. First 
bv rotating the inner grid to the specified angle of incidence and setting the oncoming 
flow to zero degrees. Second by rotating the flow to the angle of incidence and 
leaving the inner grid at zero angle relative to the outer grid. The computed pressure 
coefficients and boundary layers were the same for both cases. 

For the low Mach number viscous solutions, no flux limiting was applied. 
[t is seen that the computed results closely agree with the experimental data. At the 
higher angles of incidence the suction peak is not exactly captured. This is probably 


due to lack of grid resolution at the leading edge of the airfoil. 
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Figure 4.15: Viscous Computed Solution Compared to Experiment for a = 3.27 
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Figure 4.16: Viscous Computed Solution Compared to Experiment for a = 4.97”. 
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Figure 4.17: Viscous Computed Solution Compared to Experiment for a = 6.69°. 
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Figure 4.19: Viscous Computed Solution Compared to Experiment for a = 3.38°. 
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Figure 4.20: Viscous Computed Solution Compared to Experiment for a = 9.27°. 
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Figure 4.21: Viscous 
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Figure 4.22: Viscous Computed Solution Compared to Experiment for a = 10.99”. 
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Figure 4.23: Viscous Computed Solution Compared to Experiment for a = 11.99%. 


B. RAMP MOTION SOLUTION 

The unsteady solution for a ramp motion from a = 0? to a = 30° at Mo = 
0.3, Re = 2.1 x 109 and reduced frequency , k = 0.0127 was obtained on both a 
91 х 56 point inner grid and a 181 x 56 point inner grid. The pitch rate for the ramp 
motion Ë is defined as & = ac/2U,.. The computed lift response is compared with 


the experimental measurements of [Ref. 6] in Figure 4.24. 
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Figure 4.24: Comparison of the Measured and Computed Lift for the Ramp Motion 


Both the coarse and fine grid solutions closely predict the measured lift. How- 
ever, at the higher angles of attack, the finer grid gives higher lift. The computed 
surface pressure coefficient distributions at several angles of incidence are compared in 
Figures 4.25 - 4.38. Experimental surface pressure coefficients were available for the 
following angles of incidence: 2.94°, 5.84°, 8.91°, 11.76°, 15.5°, and they are displayed 
as diamonds in the figures. The computed surface pressure coefficient distribution 
is in good agreement with the measured data over the entire incidence range. The 
computed flowfield at the maximum angle of incidence, a = 15.5°, is mostly attached. 


A small separated flow region exists at the trailing e..ze region only. At a higher angle 
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of attack, a = 17.0%, the computed solution shows the development of the dynamic 


stall vortex in the leading edge region. 
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Figure 4.25: Computed ramp solution at a = 1.86°. 
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Figure 4.26: Computed ramp solution compared to experimental data at a = 2.94°. 
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Figure 4.27: Computed ramp solution at a = 4.14°. 
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Figure 4.28: Computed ramp solution at a = 4.87°. 
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Figure 4.29: Computed ramp solution compared to experimental data at a — 5.357. 
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Figure 4.30: Computed ramp solution at a = 6.72°. 
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Figure 4.31: Computed ramp solution at a = 8.02". 
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Figure 4.32: Computed ramp solution compared to experimental data at a = 5.91”. 
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Figure 4.33: Computed ramp solution at a = 9.85”. 
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Figure 4.34: Computed ramp solution at a = 10.50". 
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Figure 4.36: Computed ramp solution at a = 12.84”. 
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Figure 4.37: Computed ramp solution at a = 13.89°. 
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1. Boundary Layer Comparisons For Ramp Motion 

Early on in the development of the software, the computed pressure coef- 
ficients agreed quite well with the experimental data. As test cases at higher angles- 
of-attack were tested. it was discovered that the flow was not separating from the 
trailing edge of the airfoil. The probiem was discovered by investigating the boundarv 
layer profiles. À problem in the turbulence model was discovered and easily fixed. 
The lesson here is that although the computed and experimental pressure coefficients 
agree quite well it is important to check that the boundary layer profiles are reason- 
able. In figures 4.39 through 4.48 the computed boundary layers are compared to 
boundary layers computed by the interactive viscous inviscid boundary layer method 
of reference [22]. 

The comparisons are made at z/c = .5 and at r/c = .9 for angles-of- 
attack ranging from 2.94° to 15.5°. The computed boundary layer profiles compare 
quite well up to 15.5?. At r/c = .9 for the angle-of-attack of 15.5?, the comparisons 
diverge. This happens because of the different turbulence models used. At this angle- 
of attack, the flow is separating at the trailing edge so the turbulence models used 


become very important. 
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Figure 4.39: Comparison of computed boundary layer with an interactive boundarv 
laver program at the 90% chord for a = 2.94 degrees. 
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Figure 4.40: Comparison of computed boundary layer with an interactive boundary 
layer program at the 90% chord for a = 2.94 degrees. 
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Figure 4.41: Comparison of computed boundary layer with an interactive boundary 
layer program at the 50% chord for a = 5.85 degrees. 
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Figure 4.42: Comparison of computed boundary layer with an interactive boundary 
laver program at the 90% chord for a = 5.85 degrees. 
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Figure 4.43: Comparison of computed boundary layer with an interactive boundary 
layer program at the 50% chord for a = 8.91 degrees. 
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Figure 4.44: Comparison of computed boundary layer with an interactive boundary 
layer program at the 90% chord for a = 8.91 degrees. 
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Figure 4.45: Comparison of computed boundary layer with an interactive boundary 
layer program at the 50% chord for a = 11.77 degrees. 
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Figure 4.46: Comparison of computed boundary layer with an interactive boundary 
layer program at the 90% chord for a = 11.77 degrees. 
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Figure 4.47: Comparison of computed boundary layer with an interactive boundary 
layer program at the 50% chord for a = 15.55 degrees. 
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Figure 4.48: Comparison of computed boundary layer with an interactive boundary 
layer program at the 90% chord for a = 15.55 degrees. 
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2. Computed Ramp Flow Details 

In this section. flow features of the computed ramp motion are shown 
| їеэс features are demonstrated by the density contour lines. Mach contour lines. 
ortieity contour lines and mass-flux contours. The contour lines represent areas ol 
ilie flowfield, for which the parameter of interest. remains constant. Clustered contour 
lies represent areas where the parameter is rapidly changing. 

In figures, 4.49 through 4.64, density. Mach number, vorticity and mass 
Mem6ours are displayed for a = 1.86” to a = 17.50%. In figures 4.65 through 4.110 
the Alach number contours are replaced bv pressure contours. This was done becanse 
Ile strong gradients, caused by the vortices, tend to not give clear Mach number 
miormation. 

In figure 4.63, a = 16.50”. the beginning of a vortex is visible near the 
leading edge. At a = 17.00%, in figure 4.64 the vortex has moved 5% of chore 
downstream, and another vortex is starting at the leading edge. At a = 19.00 
ile first two vortices merge into one stronger vortex located near the 50% chord 
point. [wo smaller vortices are also clearly visible in the mass-flux contours of figure 
1.05. Тһе original vortex now detaches from the surface of the airfoil at 19.50%. .\ 
a = 20.50%, a second vortex is shed from the leading edge area, and a the vortex shown 
at [9.5? has moved to the trailing edge. At a = 22.50°. the trailing edge vortex has 
heen shed downstream. The whole process seems to repeat itself at a faster rate. as 
cam be seen in the remaining figures as the airfoil ramps up to 30°. The dynamic stall 
phenomenon is observed in terms of density, pressure and vorticity. [he nose-down 
pitching moment is caused by the vortex sitting at the trailing edge. It is importan 
ili the vortices do not dissipate and do not get distorted as they pass through the 


wal interface. 


л 
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Next the solution was continued to angles beyond stall. Traces of alternate 
vortex shedding from the trailing edge can be seen at a = 27°. At high angles ol 
incidence the alternating vortex shedding is very well demonstrated in figures 4.91 
through 4.110. These computed solutions are in general agreement with the findings 
of the experimental investigations of Chandrasekahara et al. Also it can be seen that 
the zonal method developed is capable of computing unsteady flows at very high 


aneles-of-attack showing at least qualitative agreement with the experiment. 
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Figure 4.52: Ramp Motion Flow Details, Ma = .3, k = .0127, Re = 2.7 x 10°, 
üt o SO. 
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Figure 4.56: Ramp Motion Flow Details; MM = eG eS 104 
a= 0 Sor. 
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Figure 4.58: Ramp Motion Flow Details, Mj, — .3, k = .0127, Re = 2.7 x 10°, 
рте. 
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Figure 4.62: Ramp Motion Flow Details, M4, — .3, k = .0127, Re = 2.7 x 10°, 
а= О 
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Figure 4.68: Ramp Motion Flow Details, M = .: => 
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Figure 4.70: Ramp Motion Flow Details, M., = .3, k = .0127, Re = 2.7 x 10°. 
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Figure 4.74: Ramp Motion Flow Details, M. 
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Figure 4.75: Ramp Motion Flow Details, Mo. = .3, k = .0127, Re = 2.7 x 10°, 
aa 2.00". 
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Figure 4.76: Ramp Motion Flow Details, Mo, = .3, k = .0127, Re = 2.7 x 10°, 
т = 25100. 
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lure 4 Ramp Motion Flow Details, Mo. = .3, k = .0127, Ве = 2.7 х 10°, 
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Figure 4.84: Ramp Motion Flow Details, Mo =. 
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mure 4-87: Ramp Motion Elow Details. WM. = .3, 6 = 0127 Re = 2.7 x 10°. 
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Inoue 4.97: Ramp Motion Flow Details. = .3, k = .0127, Re =.2.7 x 10°. 
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ENLOSCILLATORY MOTION SOLUTION 

The unsteady solution for a periodic oscillatory motion, given by a(t) = 4.86 + 
በጠ ዘጠ ን 1]. :206.)/:,ሙ= 0.6. ሸሮ. = 4.6 x 10%. with a reduced frequency of k = 0.16 
was also obtained. Here the reduced frequency is defined as KF = we/l’,,. The flow fo 


lis motion is initially purely subsonic; but. as the angle of attack increases to abou! 


if) = 5°. supersonic flow conditions are encountered at the leading edge region and а 
‘ransonic shock forms. This shock is present during the upstroke until the maximun: 
anule Of attack is reached and during the downstroke up to about a(t) = 5.0%. The 
computed and measured lift and pitching moment response are compared in Figs. 


Bernd 1.112. respectively. 


a(t)=4.86+2.44sin(wt), M=.6, k=.16, Re=4800000 
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Figure 4.111: Comparison of Measured and Computed Lift for Oscillatory Tes! 
Motion 

The computed lift and pitching moment coefficients are in close agreement wit: 
ili ıneasured values. The computed surface pressure distribution is compared with 
the measurements of reference [6] for two angles during the upstroke and two angle- 
during the downstroke in Figs. 4.113 through 4.116 The computed surface pressure 1- 


i! better agreement with the measurements at the lower angles of incidence (a = 5.95 


7 


up and a = 5.11? down). At higher incidences (Fig. 4.114 and 4.115). the agreement 
scteriorates in the region around the shock. The global view of the computed density 
lied shows that the density contours smoothly cross the zonal interface for the case 


‘here a shock exists. 


a(t)=4.86+2.44sin(wt), M=0.60, k=0.16, Re=4800000 
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Figure 4.112: Comparison of Measured and Computed Moment Coefficient for 
Oscillatorv Test Case 


a=5.95 deg. up, M=0.60, k=0.16, Re=4800000 
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Figure 4.113: Comparison of the Measured and Computed Unsteady Surface Pres 
oe. oelficient of Oscillatory Motion. a = 5.93% upstroke. 
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igure 4.114: Comparison of the Measured and Computed Unsteady Surface Pres: 
sire Coefficient of Oscillatory Motion. a = 6.97° upstroke. 
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a=6.57 deg. down, M=0.60, k=0.16, Re=4800000 
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Figure 4.115: Comparison of the Measured and Computed Unsteadv Surface Pres- 
sure Coefficient of Oscillatory Motion. a = 6.57° downstroke. 
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Figure 4.116: Comparison of the Measured and Computed Unsteady Surface Pres. 
sure Coefficient of Oscillatory Motion. a = 3.95° downstroke. 


V. CONCLUSIONS 


A solution procedure suitable for steady and unsteady compressible flow solu- 
lions using zonal overlapped grids was developed. Simple weighted averaging wa~ 
used at the overlapped zonal interfaces. Steady and unsteady. inviscid апа viscous 
(low solutions for subsonic and transonic flows over airfoils were presented to validate 
he zonal grid approach. 

The inviscid solutions presented show the overlapped zonal intertace’s ability tu 
pass flow properties without distortion. [t was found that for the inviscid test cases 
the location of the zonal interface is not important. In fact. as the zonal interface 
moved closer to the airfoil, while keeping the number of grid points constant. the 
pressure coefficient prediction actually improved. This is due to more grid poini- 
опо Clustered close to the airfoil. This test case had strong shocks on the upper ano 
«ver surface of the airfoil, and as the grid points were moved closer to the airfoil. 
‘lie computed solution tended to have bigger oscillations near the shock. 

le steady viscous test cases showed that using the Baldwin-Lomax turbulence 
model with the present approach gave accurate results for attached and mildly sep. 
arated flow over stationary airfoils. This case demonstrates one of the advantage~ 
ol the present approach, specifically, that solving the inviscid equations on the oute: 
«rid and the viscous equations on the inner grid gives good results. The flow variables 
were also passed smoothly through the zonal interface. 


the ramp case again showed good agreement with experimental data. With 


Ihis case another advantage of the present approach was displayed. The imner gri 
was rotated with the airfoil to the new angles of attack. Ihe high order accurate 


«cheme enabled the software to convect vortices, and at high angles. these vortices 


convected up to 3 chord lengths. 

The final test case was an oscillatory one. For this test case, supersonic flow was 
encountered on the upstroke at the leading edge region which produced a transonic 
shock. The computed lift and pitching moment coefficients are in close agreement 
with the measured values. At the higher angles of attack the agreement with measured 
data deteriorated in regions near the shock. 


The following are some recommendations based on this study. 


1. It has been shown that the zonal grid approach can be used to study unsteady 
viscous flows. A systematic comparison with other codes should be conducted 


in order to quantify the efficiency of the present approach. 


2. Presently, there exists no easy way to generate zonal grids. In this study the 
grids were generated separately and then refined through several other pro- 
grams. The development of a new or the modification of an existing software 
package such as GRAPE is recommended. The user should be able to specify, 
along with all the usual information, airfoil shape, the location of the overlapped 


region, the number of overlapped cells and the size of the outer grid. 


3. An advantage of the present approach is that it can be extended to multiple 
inner and outer grids. A flow solver needs to developed to take advantage of 


this so that effects of oscillating bodies in relative motion can be studied. 
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